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Abstract 

We consider a mathematical model of thermoacoustic tomography and 
other multi-wave imaging techniques with variable sound speed and atten- 
uation. We find that a Neumann series reconstruction algorithm, previ- 
ously studied under the assumption of zero attenuation, still converges if 
attenuation is sufficiently small. With complete boundary data, we show 
the inverse problem has a unique solution, and modified time reversal 
provides a stable reconstruction. We also consider partial boundary data, 
and in this case study those singularities that can be stably recovered. 



1 Introduction 

A multi-wave imaging technique couples a high-resolution imaging modality 
with a high-contrast imaging modality through a physical mechanism. In medi- 
cal imaging, ultrasound waves are often used, as they can provide sub- millimeter 
resolution in tissue. The typical examples of such techniques are thermoacous- 
tic tomography (TAT) and photoacoustic tomography (PAT). The former uses 
low-frequency microwaves as the contrast modality; the latter uses near-infrared 
light. Both are coupled to ultrasound through the photoacoustic effect. Acous- 
tic pressure is then recorded using transducers placed along the boundary of the 
region of interest. 

To reconstruct an image from the acoustic data, one solves two coupled 
inverse problems. First, one reconstructs the ultrasound source inside the region 
using the acoustic data on the boundary. Then, the ultrasound source is used 
to reconstruct the absorption of the contrast modality within the tissue. We do 
not consider the second problem here. In the case of either TAT or PAT, it is 
called cither quantitative thermoacoustic tomography (QTAT) or quantitative 
photoacoustic tomography (QPAT). 

The problem of determining the source that induces a measured acoustic 
pressure at the boundary of a region has been well-studied under the assump- 
tions of constant speed and zero attenuation. In this case, the problem is es- 
sentially the inversion of a spherical mean operator [TJ[7] . However, in tissue 
these assumptions do not hold. Variations in either acoustic speed (TO] or at- 
tenuation 4-61IT7] can cause artifacts in reconstruction. This has spurred the 
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development of reconstruction methods using models of ultrasound propagation 
that include variable acoustic speed and attenuation. 

The method of time reversal provides such a reconstruction [8,9. In this 
method, one fixes a time T (possibly infinite) and solves a mixed boundary 
problem for the acoustic wave equation with zero Cauchy data at time T and the 
measured data along the boundary for t € (0, T). This exploits the fact that the 
acoustic wave equation is invariant under time reversal, i.e., t M> — t. To ensure 
the measured data is compatible with the Cauchy data, a cut-off function is 
often used. This technique has also been extended to the thermoelastic equation 



It is also possible to choose Cauchy data that is compatible with the mea- 
sured data, thereby avoiding the use of a cut-off function. In [22] , they demon- 
strate that a modified time reversal of this type provides a stable reconstruction 
in the acoustic wave model. It has also been extended to the elastic wave equa- 
tion [25] . One goal of this work is to extend modified time reversal to a model 
including variable attuenation. An advantage of modified time reversal is that 
under some circumstances it yields an explicit Neumann series for the ultra- 
sound source. This Neumann series was studied numerically in [18] and was 
shown to be an improvement over time reversal. 

There are several commonly-used models of attenuation in biological tissue 
[HJ[T31[2D]- Generally speaking, the attenuating effect of a medium depends 
upon the frequency of the wave propagating through it, with higher frequency 
waves experiencing more attenuation. A common approximation posits a power 
law between frequency and attenuation [TBIEllEB] . Here we take the damped 
wave equation as a simple model of ultrasound attenuation. We follow the 
undamped model of [3] in assuming that the ultrasound source is pulse of the 
form f(x)6'(t). This approximation can be made in the case when the contrast 
media is an electromagnetic wave, as it is in both TAT and PAT. Writing the 
damped wave operator as D a = df + adt — c 2 A, we then consider solutions of 



in the sense of distributions. Here c is a smooth, positive function representing 
the sound speed. If f2 is our region of interest, we assume c is equal to one on 
M" \ f2. The attenuation coefficient a is assumed to be a smooth, nonnegative 
function. If we take / to be in Hq(Q), then solutions of the Cauchy problem 



are also solutions of (TjQ) when extended by zero to (-co, 0) x W 1 . To show this, 
take u to be a smooth solution of @ on K" +1 . Then we may consider the 
distribution H(t)u(t,x), where H(t) is the Heaviside function. We have that 



Urn]. 




(i) 



D a u 

u\t=Q 

d t u\ t= o 



in (0,oo) x K", 
/, 

-af. 



(2) 



D a (Hu) 



uS' + 2(d t u)8 + au8 + (D a u)H. 
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The last term is zero, as O a u = 0. Let 4> £ C ( ?°(R n+1 ) be a test function. Then 

(n a (Hu),<t>) = [ [-{d t u)^-u{d t <l>) + 2(d t u)(f> + au(p]\ t=0 dx, 

fd t (j>\t=odx, 
= (fS'^), 

which implies that u satisfies ([T]). This justifies the choice of boundary conditions 
in ©. 

Let r = dQ. be the boundary of region of interest. If we have access to 
pressure data on the entire boundary, we can model the measurements with the 
operator A a f — m|(o.t)xt- For the partial boundary case, we will take a cut-off 
function \ supported on (0, T) xL', with V C T, and model the partial boundary 
data with x^-a- As (/, — af) is in the energy space H = Hq(Q) xL 2 (Q), it follows 
from [15] that A a f is a bounded operator from H D (Q) to ^((O^T) x V). The 
inverse problem is to reconstruct / from the measured data h = A a /. In the 
damped case, the modified time reversal of [22] suggests that we construct a 
potential pseudoinverse A a h from the backward problem with mixed boundary 
data, 

D a v = 0in(0,T)xfi, 
v\t=T = 4>, / q\ 

d t v\t=T = 0, [6) 



(o,T)xr 



h. 



Here </> satisfies Acj) = on and </>|r = h\t=r- By the trace theorem, h\t=T S 
iJ 1 / 2 (r), so the problem is well-posed. We then define A a h = v\t=o- 

Our first goal is to determine whether modified time reversal is an effective 
reconstruction method for the damped wave model. In particular, we are in- 
terested in the error operator K a = A a A a — I. In the undamped case, under 
some geometric assumptions, Kq is a strict contraction |22) . This implies that 
the Neumann series 

converges. This series can be used to numerically reconstruct /, and in |18) it 
was shown to be an improvement over the usual time reversal algorithm. In the 
second section, we show that for small attenuations the above Neumann series is 
still convergent. In the third section, we consider the uniqueness and stability of 
modified time reversal in the presence of arbitrary attenuation, given complete 
data. In the fourth section, we show that in some cases partial data is sufficient 
to recover the singularities of /. In the last section we prove an a priori estimate 
for the damped wave equation with mixed boundary data, following an earlier 
estimate [15] . 
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2 Modified time reversal 



We begin with some preliminary remarks on the damped wave equation. We 
will assume that the initial data (/, —af) lies in the energy space H = Hr)(Ct) © 
L 2 (il, c~ 2 dx 2 ). Here, £fo(f2) is the completion of Cq°(Q.) under the norm 



\\f\\ 2 H D = I |V/| 2 dx. 

Jn 



By Poincare's inequality, the norm of Hu(iY) is equivalent to that of Hq(Q). 
Equip T-L with the inner product 



((a,b),(c,d)) = / Va-Vc + c 2 bddx, 



vie vctc 2 ' 

n 

which makes T-L into a Hilbert space. 

Let U C W l . For an arbitrary function w that is in both C 1 (0, T; L 2 (U)) 
and C(0, T; 7J 1 (C7)), define the local energy functional 

Eu(w,t) = - / |Vw(t,x)| 2 +c~ 2 |a t w;(t,x)| 2 dx. 

Both solutions of the wave equation and the damped wave equation are such 
functions, and it is well-known that energy is conserved in the first case and 
nonincreasing in the second case. 

For use later on, we state an a priori estimate for solutions of the nonhomo- 
geneous damped wave equation with mixed boundary data. We defer the proof 
to the last section. 

Proposition 1. Let CI be a smooth domain, and u be a solution of 



D a u 


= F 


u\ t =o 


= f, 


d t u\ t =o 


= 9 


u (o,T)xr 


= h 



subject to the compatibility condition h\t=o — f\r- Then there exists a constant 
C > such that 

sup \\{u,d t u)\\ H < Ce T IHU {||F|| L i (0 , r ;L') + ll/ll^ + MIl' + IM^} , 
te(o,T) 

provided the boundary data are sufficiently regular. 

Now, take u to be the solution of ([2]) and take v to be the corresponding 
solution of ([3]) with h = A a f. To study the error operator K a = A a A a — I, we 
study the function w — u — v, which solves the IB VP, 

D a w = 0in[0,T]xfi, 

Vj\ t =T = u\ t= T - 0, / 4 x 
dw\ t =T = d t u\ t =T, 
w\[o,T]xr = 0. 
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Recall 4> is the harmonic extension of u\{x}xr to fl, and so the boundary data 
is compatible. Let W — (w,dtw). Then ((4]) reduces to the first-order system, 



(dt + Q a )W = 0, 



where Q a is the operator 



Qa = 



-/ 
-c 2 A a 



The following energy estimate for solutions of ((4]) is used often in what follows. 
Lemma 1. Let u solve 

U a u) = Oto[0,T]xQ, 

d t u\ t =T = wi, 
, w|[o,T]xr = 0. 

Then forO<t<T, 

Ea{u,t) < e 2i - T ^ a ^E n {uj,T). 

Proof. Define the operator 

e -tQ a = e -t(Q a +||o|U/) e t||o|U/_ 

By [H Theorem X.48], if Re(Hf, /) > for all / e D(H), then e~ tH is a C°- 
semigroup of contractions for t > 0. This is true of Q a + ||a||oo as an operator 
on %. Therefore | e ^ f Qa 1 1 < e^klloo for t > 0. One can also verify directly that 



;(t) 



- p-(T-t)Qa 



(wi,a>2) is the solution of the problem in the lemma. 



Hence for < t < T, 

Enfat) < ||e-( r -*)Q (a; 1 ,w 2 )||^ < e 2T||a| 
yielding the desired energy estimate. 



□ 



This is a quantitative estimate of the nonincrease of energy in the damped 
wave equation, in a time-reversed sense. When solved backward, the damped 
wave equation gains energy at most exponentially. One consequence of this 
estimate is a bound on ||if /||iJD- 

Lemma 2. There exists a constant C > 0, independent of a, such that 

\\K a f\\ HD <(l + q|a|| 2 ) 1/2 e T||a|| -||/||H D . 
Proof. Again let W = (w,d t w). By Lemma [IJ 



\K a f\\ 2 HD < \\W(0)\\ 2 n 



< e 



2TI 



5 



Using the fact that <f> is harmonic, we see that 

\\W{T)\\ 2 H = \\(u\ t=T - ^d t u\t= T )\\ 2 u < E n (u,T). 

Hence 

||AVHL<e 2T|HU (ll/llL + IKVll| 2 ). 
By Poincare's inequality, ||c _1 a/||| 2 < C\ \a\ ^ 1 1/| \ 2 Hd , and so we have 

\\KJ\\ Hd <(l + C||a||L) 1/2 e T||Q|| -||/||H D . 

By itself, this estimate cannot show K a is a strict contraction on Hu(fl), 
this being sufficient for the convergence of the Neumann series. □ 

2.1 Convergence of the Neumann series 

Let To be the supremum of the length of geodesies in (ft, c~ 2 dx 2 ). If the metric 
is trapping, this is infinite. In [22], they show that if To < T < oo, then 
||ifo|| < 1- In that case, the Neumann series 

oo 
m=0 

converges. This series can be used as an effective numerical reconstruction 
algorithm [18]. We extend this convergence result by continuity to the damped 
wave model, given sufficiently small attenuation. 

Theorem 1. IfTo<T< oo, then there exists ao > such that for all a such 
that Halloo < ao, the Neumann series 

oo 

f = KTA a h 

771=0 

converges. 

To show this, we use the continuity of the map a n- K a . 

Proposition 2. Fix ao > 0. Then there exists C > such that for all a with 
IN loo < 

||*o/ - KJ\\ Hd < Ca (l + al) 1 l 2 e T ^\\f\\ HD 

Proof. To prove this, we return to system of coupled PDEs defining A'o and K a , 
in order to compare their difference. Let Kof = wo|t=o, where wq solves, 

n u = Q in [0,T]xl", 
«o|t=o = /, 
<9 t uo|t=o = 0, 
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n a w Q = in [0,T]xO, 

W |t=T = U \t=T — <t>0, 
d t W \t = T = d t U \t=T, 
Wo\dQx[0,T] = 0. 

As before, <po is the harmonic extension of uo\{t}x3Q- Recall K a f — w\t=o, 
where w — v — u solves ([!]), v solves and u solves @. Define u = u — u 
and w — Wq — w. Then these functions satisfy the system of PDE 

= ad t u in [0,T]xM. n , 
= 0, 
= af, 

D Q w — ad t w in [0, T] x f2, 

w\ t= T = u\t=T — {<h—<j>), 
d t w\ t =T = d t u\ t= T, 

w\dClx[0,T] = 0. 

Our interest is in Eq(w, 0), and to determine it we make use of Proposition 
[TJ Note that (</>o — 0) is harmonic. It is also equal to zero on T, and therefore 
compatible with the Cauchy data. As before Eq(u, T) is a bound on the energy 
of the Cauchy data. By Proposition [TJ 

Ea{w,Q) <c{\\ad t w\\l 1{0tT . LHQ)) + E n (u,T)}. 

Concerning the first term, by Jensen's inequality, 

ll a ^|£i (0jT:L 2 (f2)) < llal^ / \\dtw{s,-)\\ 2 L 2 {n} dxds. 

Jo 

Then by Lemma [21 we have 

\\d t w(s,-)\\h m <^K0)<(l + q|a||L)e 2T||a|| -||/Hk- 

Hence 

\\ad t w\\ L ^ T , L2m < a (l + Cal) 1 l 2 e Tao \\f\\ HD 

To estimate the second term, we use Duhamel's principle to represent u as 
the integral 

u(t, x) = / u(s;t,x)ds 
Jo 

where u(s;t, x) satisfies 

D u(s;t,x) = in [s,i]xR™, 

u(s;s,x) = af, 
dtu(s;s,x) — adtu{s,x). 

By nonincrease of energy we have 

En(u(s;-),t) < E n (u(s; •) , s) < \\a\H [\\f\\ 2 HD + E n (u, s)] 
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and as before Eq(u,s) < (1 + C , ||a||^ )e 2T H a ll 0O ||/||| /£3 . Hence by Jensen's in- 
equality again, 

Ea{u,T)< / r ^W S ;-),r)d S <C|| a ||L(l + ||a||L)e 2T||a ll-||/|||, D . 
Jo 

Combining the contributions of both terms and Theorem [T] we have 

ll*o/ - K a f\\l D < CWaWKl + ||a|| 2 00 )e 2T ll a ll~||/||^ 

which completes the estimate. Note that the constant C does not depend on 
||a||aa- □ 

Now we return to the proof of Theorem [T] 
Proof. By the previous proposition, 

\\K a f\\ < \\K f - K a f\ \ + \ \K f\\ < ||^o/||+Ca (l + a 2)i/2 e Tao|| / || HD _ 

By assumption ||ifo|| < 1- As ao — > 0, the second term tends to zero, as C does 
not depend on ao. Hence there exists ao sufficiently small such that \\K a \\ < 1 
uniformly. □ 

If 1 1 -Ka|| < 1) then the Neumann series of Theorem Q] also implies stability in 
the case of complete data. 

3 Complete data 

In this section, we use microlocal methods to prove that under some geometric 
assumptions and with complete boundary data, modified time reversal is stable. 
This follows directly from the previous section, provided the attenuation is suffi- 
ciently small. Here we show stability without this assumption. For convenience, 
we adopt the notation of [H] for unit speed geodesies. We identify T*SI\0 with 
TD, \ 0, via the Euclidean inner product. 

Notation. Let (xo,£q) € T*il \ 0. There are two null bicharacteristics of 
the damped wave equation passing over (xo,£o)- Both project to a unit speed 
geodesic (with respect to c~ 2 dx 2 ) passing through Xq. In such a situation, we 
write 7x , £ for the unit speed geodesic with 7x ,fo(0) = x o an d Tro.^o (0) = £o- 

In this section and the following one, we assume a priori that supp / is 
bounded away from T. This can always be accomplished in practice by slightly 
enlarging the region of interest, if necessary. 

3.1 Uniqueness 

First, we show uniqueness. Define d(x, T) to be the infimum of the lengths of 
curves (with respect to c~ 2 dx 2 ) starting at x and ending at T. In the undamped 
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case, [531 Theorem 2] yields uniqueness under the assumption that SI is strictly 
convex and that T > Ti(Q), where 

= supd(x,r). 

This result exploited the invariance of the wave equation under time reversal. 
To compensate for the lack of this symmetry in the damped wave model, we 
show uniqueness under the assumption T > 2Ti(Sl). 

Theorem 2. Assume that SI is strictly convex and that T > 2Ti(f2). Let 
A a / = 0. Thenf = 0. 

Proof. Let u be the solution of the forward problem ([2]). We apply Tataru's the- 
orem [33] (actually, the special case [HJ Theorem 4]). Let U be a neighborhood 
of x e R™ \ n with U n O = 0. By assumption, u = on {0} x (R™ \ SI) and 
(0, T) x L, and satisfies the damped wave equation on the exterior of (0, T) x SI. 
Therefore u = 0on (0,T) x U. By Tataru's theorem, u = on the intersection 
of the backward light cone at (T,xq) and the forward light cone at (0,xo). Re- 
peating this argument for all xq on the exterior of S7 sufficiently close to L, we 
find that u = 0on 

(t,x):d(x,T) < |- ~-t 

Since T > 2Ti(Sl), in particular we see that w = on a small neighborhood of 
{T/2} x ft. By Proposition [U u = on (0, T/2) x SI, and so / = 0. □ 



3.2 Stability 

To show the reconstruction is stable even when attenuation is large, we turn to 
microlocal analysis. When a = 0, stability was shown in [351 Theorem 3] under 
the assumption that T was sufficiently large so that for every (x, £) 6 S*fl, 
either 7 Xi j or 72, hit T before time T. Recall from the previous section that 
To is the supremum of the lengths of geodesies in (SI, c~ 2 dx 2 ). Here we present a 
stability estimate in the special case where c -2 dx 2 is non-trapping, i.e., To < 00. 

That this assumption is not necessary for stability follows from the result 
for partial data in the next section. However, the proof in this case is simpler. 
Notice that in the non-trapping case, if T > Tq(CI), then T > 2Ti(Sl) as well, 
and so Theorem [2] implies that A Q is injective. 

In the following theorem, let x G C°° x T) be a fixed cutoff function equal 
to one on [0, T ] x L and zero in a neighborhood of {T} x L. 

Theorem 3. Assume To < T < 00. Then the following are true: 

1. A a A a : H D (Cl) -)• -ffo(Sl) is Fredholm. 

2. A a xAa = I + R a , where R a is a smoothing operator. 

3. There exists a constant C > such that 

\\f\\H D <C\\A a f\\ m . 
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Proof. First, we show A a A a is Fredholm. In the second section, we considered 
the error operator K a = f — A a A a f. It remains to show that K a is compact. 

Let u be the solution of the forward problem @. By propagation of sin- 
gularities, WF(u) is the union of null bicharacteristics lying over WF(/). Let 
h = A a f. As we assume T > To (17), WF(w) is empty over a neighborhood 
of {T} x 17. Therefore the operator / t-» (u\t = r — <fi,dtu\t=T) maps £/d(17) to 
C°°(17) x C°°(17), and hence is compact. Finally, Proposition Q] implies that the 
operator (u\t=T — 4>i dtu\t=r) >-> w\t=o is bounded from H — > Hd(Q). Therefore 
K a is compact, and A a A a is Fredholm. 

Let v be the solution of the corresponding backward problem ([3]) with h = 
xA a f, and w = u — v as before. Evidently v\ t= T = 0. Therefore w solves 
((4]) with smooth boundary data, compatible to arbitrary order. If we dchnc 
R a f = w\ t= o, then R a maps H D (Q) to C c °°(17), and A aX A a = I + R a . 

To show the first inequality, the previous equation shows 

\\f\\ HD < \\A aX A a f\\ HD + \\R a f\\ L 2. 

By Lemma [2J A a is bounded from -ff 1 (r) — > ij£>(17), and so we have 

\\f\\H D <C'(\\A a f\\ H1 + \\f\\ L i). 

To obtain the stability estimate in the theorem, we recall that since T > 2Ti(f2), 
A a is injective. So by [24j Proposition V.3.1], there is another constant C > 
such that 

||/||h d <c||A n /|| ff i. 

Therefore the reconstruction is stable. 

If in addition 17 is strictly convex (with respect to the Euclidean metric), 
the results of the next section show that singularities that propagate through 
the boundary may be stably recovered; the larger T is, the more singularities 
become visible. □ 

Concerning the second part of the theorem, in principle one could choose a 
different modification of the Cauchy data in ([3]) that would be compatible with 
h to at least first order. In that case, K a would be smoothing of at least first 
order, and no cut-off function would be necessary. 

4 Partial data 

In applications it is often only possible to obtain data on a proper subset of 
the boundary. Let f C T be relatively open. We model partial measurements 
by the operator xA Q , where x S C°°(IR x V) is a cut-off function positive on 
(0, T) x T' and zero elsewhere. In this section, we only show that "visible" 
singularities (in the sense of [21E2]) ma y be stably recovered from xA tt /, by 
proving an estimate of the form 

\\f\\H D <C(\\A a f\\ H1 + \\f\\ L2 ). 
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If in addition, we know that A a is injective, we could prove stability as in the 
complete data case. In the undamped case, [22l Proposition 2] shows uniqueness 
provided that d(x,T') < T. However, their result relies on being able to extend 
solutions of the wave equation to (— T, 0) x f2 by time reversal. It is not clear 
how to modify their proof to account for attenuation. Therefore, we constrain 
ourselves to showing stable recovery of visible singularities. 

Let K. be an open set strictly contained in f2. By analogy with the complete 
data case, let T 2 (IC,T') be the least time necessary for all G S*K. to have 

the property that at least one of y X £ and 7 x ,-$ reaches T' before time T 2 (IC, V). 
In other words, we can say that K, is visible from V. We will assume that 
T > Ta(/C, T'). Our goal in this section is to show the following: 

Theorem 4. Assume fl is strictly convex (with respect to the Euclidean metric), 
and that T2(K.,V) < T < oo. If in addition, supp/ C IC, then the following are 
true: 

1. A a xA a : H D (IC) ->• H D (Q) is Fredholm. 

2. A a \A a : _ffo(f2) Hrj{Q) is an elliptic pseudodifferential operator of 
order zero. 

3. There exists C > independent of a such that 

\\f\\ HD <Ce T ^{\\K a f\\ H , + \\f\\ L ,). 

Proof. The first part follows directly from the second part of Theorem [31 though 
one may need to increase T slightly to ensure that x = near {t — T}. 

To show the second part, we construct a geometric optics solution u of ([2]), 
and then microlocally construct a parametrix of (J3J given h — x^-a- We defer 
the proof of this to the following two lemmas. 

Finally, provided the previous part, by elliptic regularity we have the esti- 
mate 

\\f\\H D <C(\\A aX Kf\\H D + \\f\\ L i)- 

By Proposition [I \\M\mp)-+H D <fi) < Ce T ^°°. □ 

The principal symbol of D a factors as (r + c(x)\r]\)(—T + c(x)\rj\). Our geo- 
metric optics solutions will then depend on the existence of two phase functions 
solving the corresponding eikonal equations. To simplify the proof of the next 
two lemmas, we will assume that these equations are solvable on the interval 
[0,T], i.e., that no caustics form. Afterward, we will remove this assumption. 

Lemma 3. A Q : Hrj(fl) — > H (T) is the sum of two Fourier integral operators, 
A+ and A~ , with canonical relations 

{(y,r),t,x,T,£) : t = T y .± v ,x = y y ,±r,(t),T = -\j y ,± v (t)\,£= -7r(7j, )±J7 (t))}, 

where ir : T*W l — > T*T is the tangential projection. 
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Proof. We will search for a parametrix u for @ of the form 

u(t,y) = (2tt)- b ^ / e**'^A"(t,y,ri)f( V )dri, 

where <jr^ are phase functions and A 7 are order zero amplitudes with asymptotic 
expansion 

j>0 

where A° (t,y,r]) is homogeneous of degree — j in 77. For u to be a parametrix, 
we must construct (jf and A CT so that O a u € C°°, u| t= o = /, and <9tu| f= o = —a/. 
To satisfy the first requirement, we calculate □„«: 

□„u = (2tt)- w ^ / e**7[Ja +h + Jo] dq, 

where 

h = 2i[(d t r)(dtA a ) - c 2 \7 y r ■ V y A°] + iA°U a <f, 

h = u a A a . 

For the residual term to be smooth, the integrand must decay rapidly with 
respect to 77. If we impose the eikonal equation 

f T a^± = c |v^±| 
\ 0^=0 = y-v 

and assume for the moment that a solution exists for t G [0, T], then evidently 
h=0. 

We recognize the first term of I\ as the transport operator X a applied to 
A", where 

X° = 2(d t r)dt - 2c 2 \7 y r ■ V y 

To control the decay of 1% + 1% in r], we write A a as its asymptotic expansion 
and consider the coefficients separately. On Aq, we can impose the transport 
equation 

X a Al + AZU a ^ = 0. 
The subsequent terms can be taken to satisfy the lower order transport equations 

X a A a j + A^U a (jf = -□ AJ_ 1 ,j > 1. 

These equations reduce to ordinary differential equations along character- 
istics - the characteristics again being unit speed geodesies - whenever the 
eikonal equations are solvable, provided we establish Cauchy data at t = 0. To 
determine the proper Cauchy data to impose, we return to the Cauchy data of 
©■ 
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The condition u\t=Q = / implies that when t = 0, 

f(y) = (2n)- n J e^f( V )[A + +A-]\t=od V , 
and accordingly at t = 0, 

A + + A~ = l. 
The condition dtu\ t =o = — a/ implies that when t = 0, 

/ e iv ' n f(v)(-a(y))dr) = J e^/fo) [ic|?7|(-^ + + A") + <9 f (A+ + A - )] dr?. 

Therefore at i = 0, 

ic|?7|(-A + + A~) + d t {A+ + A~) = -a 

These two equations yield the following system of boundary conditions at 
t = 0: 

A-t + A) =1, ( Aq — Aq = 0, 

At + A^ = 0, { A+-A^ = -a-d t (A+ +A ), 

A+ + AJ = 0, j > 2 { Af - Af = -dtiAf^ + Af_J, j > 2. 

Note that the resulting system of equations may be solved recursively, yield- 
ing boundary data for each coefficient of the asymptotic expansion of A^. In 
particular, A^(0,y,r]) — A^(0,y,7]) = |. With this set of boundary data, 
the recurring system of transport equations stated above may be solved. This 
completes the construction of u. 

The restriction of this parametrix to (0, T) x T yields a global representation 
of A a as the sum of two oscillating integrals, up to a smooth residual error. For 
(t,x) G (0,T) x T, we have 

[Atf](t,x) = (2n)- n J e^ ± ^A ± (t,x, V )f(r i )dr ! . 

From this we can calculate the canonical relation of each Fourier integral oper- 
ator. The characteristic manifold is 

S ± = {y = d 7l <p ± (t,x,r ] )}. 

It is well-known that y = 9 7) </> ± (t, x, rf) holds when ^^(t) = x. Therefore, the 
canonical relation of A„ , restricted to T*Q xT*((0, T) xT), is given locally by the 
graph of the diffeomorphism mapping (y 0) rj ) to (T y(u ± Va ,~f y(u ± Va (T y(u±m )). □ 

Next we perform a similar construction for the backward problem ([3]), in 
order to analyze A a \A a . 

Lemma 4. A a \^-a '■ Hu(Q) — > Hd(£1) is a pseudodifferential operator of order 
zero whose principal symbol is 

O"o(y,»7) = \[x{Ty,r,,ly t 7,{ry, v )) + X(r y ,- V ,7y,- V (r y ,- V ))] 
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Proof. As in the previous lemma, we begin by constructing a parametrix v of ([3]), 
assuming again that the eikonal equations are solvable on (0, T). To proceed, we 
microlocalize the parametrix near the point (to, xo , To, £o) £ T*((0, T) x T). We 
take p 6 C£°(R n+1 ) to be a cutoff function supported in a small neighborhood 
U of (to,xo) such that p(io,^o) = 1- If T o < 0, we will take h = px^-t f an d 
assume that WF(ft) is contained in a small conic neighborhood of the point. If 
r > 0, replace A+ with A~ and proceed similarly As Cauchy data for ([3]), we 
will take v\ t= r — d t v\ t= T = 0. The boundary data is smoothly compatible, as 
we assumed that \ is supported away from {t = T}. 

Under these considerations, we search for v solving ([3]) up to smooth error 
with the above initial data. We assume v is of the form 



as suggested by the forward problem. 

Propagation of singularities implies that the wavefront set of v should be 
near the two null bicharacteristics whose projection onto Tt Xo -,((0,T) x V) is 
( t Oj£o)- The assumption of strict convexity implies that there are two distinct 
geodesies whose tangent at xq projects to £o! one pointing inward, and one 
pointing outward. 

Let (yo,Vo) be related to (to,xo,TQ,£o) under the canonical relation of A+. 
Therefore, the unit speed geodesic "f yo ,ri is such that xo — 7»0)'7o( 7 Vo>»7o 

). Evi- 
dently, this is the geodesic pointing outward. 

The other geodesic is the reflection of this one along T, according to the laws 
of geometric optics. Singularities near it propagate further, perhaps reflecting 
off the boundary at most finitely many more times (as the boundary is strictly 
convex), until they intersect {t = T}. Since we specified that the Cauchy data 
is zero there, these reflected broken geodesies carry no singularities. 

Since only the null bicharacteristic corresponding to jyg^o carries the singu- 
larities of h, we need only construct ip and B in a small conic neighborhood of 
this bicharacteristic. The boundary data is of the form 



This holds in particular if we impose — dtv = c(x)\ y 7 x ip\. Then if) and <f> + are 
given by the method of characteristics in a small conic neighborhood of the 
bicharacteristic. They have the same characteristics and the same boundary 
data, so in this neighborhood they are equal. 

Similarly, B solves the same recurring system of transport equations as A + , 
though with different initial data; on U, B — p\A + . Since the first transport 





For v to be a parametrix, we must have 



ip\u 



c 2 (z)|V x Vf 
<P + \u- 
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equation is homogeneous, we have, 

B (0,y ,r]o) = p(t ,x )x(to,xo)A£ (0,y ,r} ) = ^x(Ty , vo ,ly ,r, {T yo ,r, )). 

We repeat the construction for points in T*((0, T) x V) with To < 0, as men- 
tioned before. Passing to a microlocal partition of unity, the above constructs 
v as a parametrix for ([3]) with boundary data h — x^-af and zero Cauchy data 
at {t = T}. Restricting v to {t — 0} yields a representation of A a \A a as a 
pseudodifferential operator, whose principal symbol is Bq. After accounting for 
both contributions to Bq(0, t]q) from singularities propagating along 7j, 0!?7o 
and 7, /0! _, )0 , we see that the principal symbol is 

tVOtVo) — „ [X^yotVo > 7yo,»?o ( T yo,Vo )) +x(Ty ,- no ,-fy ft ,-no( T yo,-Vo))}- 
Therefore A a xA. a is of order zero. □ 

Remark. The previous two lemmas assumed that the eikonal equation was solv- 
able on (0, T) x fi. In general, solutions with Cauchy data on {t = 0} exist 
only on a small interval [0, £i]. To continue past {t — t\}, we use the previ- 
ous parametrix to obtain Cauchy data on this hyperplane. Then, we solve the 
eikonal equations with boundary data cj)f\t=t 1 = y ■ T), and the transport equa- 
tions with Cauchy data given by the previous parametrix. This will be solvable 
on some interval (t%, £2]- By compactness, there exists a positive lower bound on 
the length of the interval on which the eikonal equation is solvable. Therefore, 
after repeating the construction at most finitely many times (when T < 00), 
we obtain a parametrix on (0,T) x Q. The backward parametrix is then also 
constructed in a similar way on the same intervals. 

We now conclude the proof of Theorem 01 According to Lemma HJ A a x^a 
is a pseudodifferential operator. Its symbol is of order zero by construction. If 
K, is visible from V , then the principal symbol comprises two terms; both are 
nonnegative, and at least one is positive. Hence A a \^ a is elliptic. If JC is not 
visible, then A a \A a is still a pseudodifferential operator, but it is only elliptic 
at those points of T*fl \ which are visible from the boundary. 

Finally, we note that even if K, is not visible from V , A a x^ a remains a pseu- 
dodifferential operator, and is still Fredholm by the considerations of Theorem 

El 

5 A priori estimate 

We return now to the proof of PropositionQ] We follow [KJ Theorem 4.2], using 
the method of multipliers to develop an a priori estimate. There, they studied 
operators of the form &f — A, with A uniformly elliptic and in divergence form. 
For our purposes here we modify their proof to account for the lower order term 
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in the damped wave equation. Recall u is a solution of 

D a u = Fin(0,T)xtt, 

u\t=o = /, 

d t u\ t =o = g 

u l(o,T)xr = h 

subject to the compatibility condition h\t=o = f\r- 

Proof. The proof has three steps. First, we establish an a priori estimate of 
the form above, except that it also depends on the L 1 (0, T; L 2 (r))-norm of 
the normal derivative of u, along with and ||9tu||i2. Next, we remove 

the first dependency by estimating the normal derivative. At this point the 
estimate contains a term depending on the L 1 (0, T; L 2 (T))-norm of Vit, which 
we separate into its normal and tangential components. The final dependencies 
are removed via Gronwall's inequality. 

Let n be the outward normal vector field of T. As T is smooth, we can extend 
n to a smooth vector field in a neighborhood of T. Using a cut-off function, we 
can then extend it to a smooth vector field on a neighborhood of f2, which we 
will also write as n. The result does not depend on the particular extension 
chosen. 

Step 1: We begin by fixing s G (0,T) in the PDE and multiplying through 
by 2c~~ 2 dtu, integrating with respect to fl. 

2c- 2 (dfu)d t u - 2{Au)u t + 2a C - 2 (d t u) 2 dx = / 2 C - 2 (d t u) 2 F dx. 

We apply the identity 

2c- 2 {d 2 u)d t u - 2(AnH = j t (Hc-^^H 2 , + ||Vu|||,) - 2V • (d t u)(Vu) 

to the first two terms. Rearranging, we have 

4- (\\ c ~ l9 t u \\% + \N u \\h) = / 2V-{d t u)(Vu)-2ac- 2 {dtu) 2 + 2c-~ 2 {d t u)Fdx. 

The second term can be simplified with the divergence theorem. After integrat- 
ing from to s and freely applying Young's inequality, we have 



Wc-'dtu^Wh + l|Vn( S )||| 2 < C |||F||| 1(0iT . i2) + Il/H^ + \\g\\l* + \\h\ 



2 



\Vu-n\ 2 dSdt 

r 



(1 + Halloa) / ||c 1 {d t u)\\ 2 L 2 m dt 



i 2 (n) ( 

This agrees with the analogous estimate [15l 4.12]. Here C is a generic constant 
that is independent of a. 
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Step 2: To establish the needed estimate on the normal derivative, we repeat 
the process above, multiplying through by c~ 2 (Vu • n) and integrating over 
(0,s)xfL 



f ( [cr 2 d 2 u- Au + c~ 2 ad t u}(Vu ■ n)dxdt = f f c7 2 F(Vu ■ n) dx 
Jo Jn Jo Jn 



dt (5) 



The first two terms reduce nearly the same as the analogous estimates [TH 4.13] 
and [13 4.16]: 

C - 2 (d?u){Vu-n)dxdt= [ cT 2 (<9 t u(s, -))(du{s, •) • n) - g{V f ■ n) dx 
u -fa Jn 

[ J(V • {c- 2 n))(d t u) 2 dxdt 
o Jn 2 

dt{c~ 2 ){dtu){'S/u ■ n) dxdt 
(V(c- 2 )-n)(d t h) 2 dSdt, 



o Jn 

1 

o Jr 2 



/ (Au){Vu ■ n) dx dt = [ ( \\V u\ 2 - \V u ■ n\ 2 dS dt 
Jo Jo Jr 2 



+ / / <(Vn)(Vu), Vu)dxdt. 
Jo Jn 



Here we use the fact that 



i|Vu| 2 - |Vtt • n| 2 = -\Vu ■ n^\ 2 - -|V« • n| 2 < -|Vit • n\ 2 
Isolating this term to the left-hand side, we obtain, 



|Vu • n\ 2 dSdt < / c~ z (dtu(s,x))(Vu(s,x) ■ n) -g(Vf-n)dx 
ii .'r Jn 

■ c- 2 n)(d t u) 2 - (d t c- 2 )(d t u)(Vu-n)dxdt 

-2\ „\/o i,\2 



o Jn 

1 



o 



-(V(0 ■ n){d t hy dS dt 



+ ((Vn)(Vit), Vu) (it 

Jo in 



+ / / c a(d t u)(Vu ■ n) dx dt 

Jo Jn 

- [ [ cr 2 F(Wu-n)dxdt. 

Jo Jn 



On the right hand side, we majorize with the absolute value and apply Young's 
inequality to several terms. This yields the desired estimate for the normal 
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derivative in terms of the boundary data and the I/ 1 (0, T; L 2 (f2))-norms of 
c~ 1 {dtu) and |Vu|. Notice it is only these last two terms that depend on | |a| |oo , 
and that the dependence is at most linear. 

Step 3: We begin from the inequality derived from applying the estimate 
for the normal derivative to the estimate from step 1. 

\\{u,d t u)(s)\\ 2 H < c{||F||| 1(0iT . £a) + H/ll^ + \\g\\l* + \\h\\ 2 Hl } 

+ (l + IM|oo) f \\{u,d t u){t)\\ 2 H dt 
Jo 

By Gronwall's inequality, for all < s < T we have 

\\(u,d t u)(s)\\ 2 n < Ce 2s IHU \\1 1{0T . L2) + \\f\f Hl + M % + \\h\\h}- 

Finally, we take the supremum over s € [0, T] to obtain the estimate of Propo- 
sition [TJ While we assumed that Q,c, and a were all smooth, the proof shows 
that their being C 1 suffices. □ 
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